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Abstract This paper provides an overview and a protocol of molecular clock dating using 
MrBayes. Two modern approaches, total-evidence dating and node dating, are demonstrated using 
a truncated dataset of Hymenoptera with molecular seguences and morphological characters. The 
similarity and difference of the two methods are compared and discussed. Besides, a non-clock 
analysis is performed on the same dataset to compare with the molecular clock dating analyses. 
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1 Introduction 


MrBayes is a software for Bayesian phylogenetic inference (Huelsenbeck and Ronquist, 
2001; Ronguist and Huelsenbeck, 2003) and has-become widely used by biologists (Van 
Noorden et al., 2014). Many new features have beenamplemented since version 3.2 (Ronquist 
et al., 2012b), including species tree inference under the 4nulti-species coalescent model (Liu, 
2008; Liu et al., 2009); compound Dirichlet priors for branch lengths'(Rannala et al., 2012; 
Zhang et al., 2012); marginal model likelihood estimation using stepping-stone sampling (Xie 
et al., 2011); topology convergence diagnostics using the average standard deviation of split 
frequencies (Lakner et al., 2008); BEAGLE library support (Ayres et al., 2012) and parallel 
computing using MPI (Altekar et al., 2004). 

In addition to these features, the focus of this study is its new functionality of divergence 
time estimation using node dating (Yang and Rannala, 2006; Ho and Phillips, 2009) and 
total-evidence dating (Ronquist et al., 2012a; Zhang et al., 2016) approaches under relaxed 
molecular clock models (Huelsenbeck et al., 2000; Thorne and Kishino, 2002; Drummond 
et al., 2006; Lepage et al., 2007). Both dating techniques are implemented in the Bayesian 
framework, such that the diversification process and the fossil information are incorporated 
in the priors. However, they do have significant differences. In node dating, only molecular 


sequences from extant taxa are used, and one or more internal nodes of the extant taxon tree 
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are calibrated by user-specified prior distributions, typically derived from the fossil record, 
to estimate the ages of all other nodes in the tree. In total-evidence dating, extant and extinct 
taxa are all included in the tree, and morphological characters are coded for fossil and extant 
taxa in a combined matrix with molecular sequences. The age of each fossil is assigned a prior 
distribution directly. 

Several steps are involved in a Bayesian molecular clock dating analysis, importantly 
including partitioning the data, specifying evolutionary models, calibrating internal nodes 
or fossils, and setting priors for the tree and the other parameters. In this paper, I introduce a 
general clock dating protocol that could potentially be applied to a wide range of taxonomic 
groups, using a truncated Hymenoptera data as an example. 

The original data analyzed by Ronquist et al. (2012a) and Zhang et al. (2016) includes 
60 extant and 45 fossil Hymenoptera (wasps, ants, and bees) and eight outgroup taxa. The 
data were divided into eight partitions as follows: 1) morphology, 2) 12S and 16S, 3) 18S, 4) 
28S, 5) I "and 2codon positions of CO1, 6) 3" codon positions of COL (but not used in the 
analyses), 7) 1“ amd 2#eodon positions of both copies of EFla, and 8) 3" codon positions of 
both copies of EFla. These are based on morphological characters, mitochondrial and nuclear 
genes. The full data takes days tố rum. For illustrative purposes, the data is truncated into ten 
extant taxa (nine Hymenoptera and one Raphidioptera) and ten fossils, with 200 morphological 
characters, 100 sites from 16S, and 210 sites from EFla. Below, the dataset is analyzed by 
both total-evidence dating and node dating approaches,under diversified sampling of extant 
taxa (Höhna et al., 2011; Zhang et al., 2016), followed bysa.non-clock analysis under the 
compound Drichlet prior for branch lengths (Rannala et al., 2012; Zhang et al., 2012). 


2 Protocol 


2.1 MrBayes installation 


The program MrBayes is available from GitHub (https://github.com/NBISweden/ 
MrBayes), including pre-compiled executables for macOS and Windows, and source code for 
all platforms. The current version is 3.2.7, which is used here. The truncated dataset hym. nex 
is in the doc/tutorial folder within the release. For convenience, it is recommended 
to put it in the same directory as the executable (e.g., named mb in macOS/Linux or mb. exe 
in Windows). The file hym. nex in the NEXUS format can be opened by a text editor or 
a matrix editor such as Mesguite (https://Www.mesguiteproject.org). The data matrix is at 
the beginning, including morphological characters and molecular seguences. Following the 
data block, users can write MrBayes commands within the mrbayes block (each ends with a 
semicolon). These commands will be executed automatically when the data is read in. The 
texts within a pair of sguare brackets are comments and will be ignored by the program. 

In terminal (macOS/Linux) or command prompt (Windows), navigate to the folder 
containing the executable and data file using the cd command, and launch MrBayes using . / 
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mb (macOS/Linux) or mb .exe (Windows). The following header will appear. 
MrBayes 3.2.7 x86 64 

(Bayesian Analysis of Phylogeny) 

Distributed under the GNU General Public License 


MrBayes > 


The prompt at the bottom means that MrBayes is running and ready for your commands. 
Simply use 
execute hym.nex 


to read in the data. 


2.2 Data partitions 


When the data is read in, the commands for defining the partitions are also executed. 
chabset MV = 1-200 
chege ES = 201-300 
charset W M = 301-510 
charset Eflal27=m301-510\3 302-510\3 
charset Efla3 ="303-51013 
partition four = 4: MW, @6sS,, Eflal2, Efla3 


set partition = four 


Here we define four partitions: the morphology (MV), 16S, 1“ and 2" codon positions of 
EF la (Eflal2), and 3" codon positions of EFla (Efta3). The name “four” is user defined and 
can be any other string. Note the semicolon (;) at the end of each command is neglected as the 
command is supposed to be typed after the MrBayes prompt (a.semicolon is necessary when 
the command is written in a file, the same below). 


2.3 Substitution models 


For the morphological partition, the Mkv model (Lewis, 2001) is used with variable 
ascertainment bias (only variable characters scored), equal state frequencies and gamma ‘rate 
variation across characters. 

lset applyto = (1) coding = variable rates = gamma 
Constant characters will be excluded automatically by the program. To be specific, use the 
following command to exclude them. 

exclude 7 31 61 83 107 121 122 133 182 183 198 
If instantaneous change is only allowed between adjacent states (e.g., 0-51 and 1+>2 but not 
0<>2), these characters are specified using this command. 

ctype ordered: 20 23 27 30 36 41 42 44 46 48 59 65 75 78 79 89 
99 112 117 134 146 157 159 171 185 191 193 196 
The other characters are thus unordered which can change instantaneously from one state to 
another. Dollo or irreversible character is not supported in MrBayes at the moment. 

For the molecular partitions, the general time-reversible model is used with gamma rate 


201906.00026v1 


chinaXiv 


Chinaxiv21EERTI 


4 


variation across sites (GTR+T) (Yang, 1994a,b). 

lset applyto = (2,3,4) nst = 6 rates = gamma 
The prior for the gamma shape parameter is exponential(1.0), which can be changed using 
prset shapepr. We keep the default here. 

Different partitions are assumed to have independent substitution parameters, thus we 
unlink them. The partition-specific rate multipliers are used to account for rate variation across 
partitions with average to 1.0. 

unlink statefreg = (all) revmat = (all) shape = (all) 


prset applyto = (all) ratepr = variable 


2.4 Relaxed clock models 


The relaxed clock models account for evolutionary rate variation over time and among 
branchéSpand it is now a standard practice to accommodate such variation in dating analyses. 
There are three relaxed clock models implemented in MrBayes: compound Poisson process 
(CPP; Huelsenbeck et al., 2000), autocorrelated lognormal (TK02; Thorne and Kishino, 
2002) and independent gammy rate (IGR; Lepage et al., 2007). However, the CPP model is 
computationally not compatible with total-evidence dating, thus we only focus on the IGR and 
TK02 models. The outcomes of using these two models are discussed below. 

The mean clock rate (mean substitution rate per site per Myr) is assigned a lognormal(-7, 
0.6) prior, with mean 0.001 and standard deviation 0.0007. 

prset clockratepr = lognorm(-7,0.6) 
The prior is chosen by comparing the age of the oldest insect fossil.with the root age estimation 
from uncalibrated clock analysis as discussed in Ronquist et al.(2012a). One might need 
to specify a different suitable distribution. There are several options fer the clock rate prior, 
including fixed, normal (truncated at zero), lognormal, and gamma. The probability 


density functions (all with mean 0900) and 
8 - — lognormal(-7, 0.6) standard deviation 0.0007) are shown in Fig, IX 
= — gamma(2, 2000) . i A i 
— normal(0.001, 0.0007) The differences in this case are minor. 
ae The relative clock rates, which are mul- 
KZ tiplied by the mean clock rate, vary differently 
= g | along the branches of the tree in different models. 
y The TK02 model assumes that the relative rate 
5 changes along the branches as Brownian motion 
sd on the log scale, starting from 1.0 (0.0 on the log 
scale) at the root. The rate at the end of a branch 
AAA Ï$ lognormally distributed with mean equal to the 
ngài o s 0904 000 rate at the beginning of the branch. Thus the rates 
Fig. 1 Probability density functions at different branches are autocorrelated. 
of lognormal, gamma, and normal distributions, prset clockvarpr = tk02 


all with mean 0.001 and standard deviation 0.0007 prset tk02varpr = exp(1) 
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The IGR model assumes that the relative rate at each branch is an independent gamma 
distribution with mean 1.0. The variance is proportional to the branch length. Thus the rates are 
independent of each other (but not identically distributed). 

prset clockvarpr = igr 

prset igrvarpr = exp(10) 
Note that when the relative rates are all fixed to 1.0 in the TK02 or IGR model, it becomes the 
strict clock model. 

Before the total-evidence dating and node dating analyses, we first define the fossil 

taxa and some constraints for later use. Note these constraints are not enforced until we set 
topologypr explicitly (see below). 


taxset fossils = Asioxyela Nigrimonticola Xyelotoma Undatoma 


Dahuratoma Cleistogaster Ghilarella Mesorussus Prosyntexis 


Psemdoxyelocerus 
constraint+HymenFossil = 2-. 
constrainteHŷmenoptera = 2-10 
constraint HoPometabola = 1-10 
constraint Tenthredińidae = 3-5 


constraint CepSirOruApo/ =) 7720 
The numbers here refer to the taxon numbersnordered as in the data matrix, while a dot (.) 
means the last taxon. 


2.5 Total-evidence dating 


In the following, we assign priors for the fossils from the geological times. This is a 
typical step in total-evidence dating, where we calibrate the fossil taxa instead of the internal 
nodes of the tree. Fossil age uncertainties are represented as uniform distributions, while fixed 
values can be used when the geological age is very certain. 

calibrate 
Asioxyela = unif (228,242) 
Nigrimonticola = unif (152,163) 
Xyelotoma = unif (152,163) 
Undatoma = unif (145,152) 
Dahuratoma = fixed (134) 
Cleistogaster = unif(168,191) 
Ghilarella = unif (113,125) 


Mesorussus = unif (94,100) 


Prosyntexis = unif(80,86) 


Pseudoxyelocerus = fixed (182) 


prset nodeagepr = calibrated 


To enable the calibrations, the last nodeagepr command must be present. 
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The speciation, extinction, fossilization, and sampling process is explicitly modeled using 
the fossilized birth-death (FBD) process (Stadler, 2010; Heath et al., 2014; Gavryushkina et al., 
2014; Zhang et al., 2016). 

prset brlenspr = clock:fossilization 


Diversified sampling (Zhang et al., 2016) is arguably suitable for such higher-level taxa, which 
assumes exactly one representative extant taxa per clade descending from time x,,, is sampled, 
and the fossils are sampled with a non-zero rate before x,„„ and zero after. A fossil can be either 
a tip or a sampled ancestor (Fig. 2). 

prset samplestrat = diversity 
By default, the sampling strategy for extant taxa is uniformly at random (Samplestrat = 
random) which is not used here, but see Zhang et al. (2016) for the application. 


A complete tree B diversified sampling 


t 


mrca 


Fig.2 The fossilized birth-death (FBD) process and diversified sampling of extant taxa 
Exactly one representative taxa per clade descending from time x,,, is sampled (blue dots). The fossils are 


sampled with a constant rate between fe and x,,, (red dots) 


mrca 


The model has four parameters: speciation rate A, extinction rate u, fossil discover rate y, 
and extant sample proportion p. For inference, we re-parametrize the parameters as'd 34 u, 
r=ul),ands=w/(u+ y), so that the latter two parameters range from 0 to 1. p is fixed to 
0.0001, based on the living number of Hymenoptera at about 100000 (10/100000 = 0.0001). 
prset sampleprob - 0.0001 


prset speciationpr = exp(10) 

prset extinctionpr = beta(1,1) 

prset fossilizationpr = beta(1,1) 
prset treeagepr - offsetexp(300, 390) 


The FBD prior is conditioned on the root age (/„„„). Here we use an offset exponential 


distribution with minimal age 300 Ma (according to the oldest neopteran fossil) and mean age 
390 Ma (the oldest insect fossil). Users should specify a reasonable root age prior according to 
the available fossil information. Several available probability densities all with mean 390 and 
minimal 300 are shown in Fig. 3. The exponential prior is more diffused than the rest three. 

An alternative tree prior that can be used in the total-evidence dating approach is the 
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uniform prior (Clock: uniform) (Ronquist 


0.012 


— offsetlognormal(300, 4, 1) 

— offsetgamma(300, 2, 0.022) 
— offsetexp(300, 0.01111) 

— truncatednormal(300, 390, 60) 


et al., 2012a). It assumes that the internal nodes 


are draw from uniform distributions and the 


0.01 


fossils are only tips of the tree (so-called tip 


> 8 
dating). Similarly to the FBD prior, the uniform 2 5 
prior is also conditioned on the root age and > E 
requires setting treeagepr. 8 z 
In order to root the tree properly, we enforce the 3 
constraint HymenFossil defined above. This Š 
forces the Hymenoptera with fossils to form a s 
paonophyletic group. ° 300 350 400 450 500 550 600 
erset topologypr = age 
CoństPai nt (HymenFossil) Fig. 3 Probability density functions of offset 
Fot the Markov chain Monte Carlo lognormal, offset gamma, offset exponential, 
(MCMC) (Metropolis et al., 1953; Hastings, and irincated normal distributions, 


8 all with mean 390 and minimum 300 
1970), we use two independent runs and four 


chains (one cold and three hot) perrun for 500,000 iterations and sample every 100 iterations. 
mcmcp nrun = 2 nchain #44 ngen = 500000 samplefr = 100 
mcmcp filename = hym.te printf», = 1000 diagnfr = 5000 
mcmc 
The output file names all start with hym. te (represented ash ym. te. *). The chain states are 
printed to screen every 1000 iterations, and the convergence.diaenostics are printed every 5000 
iterations. The mcmc command executes the MCMC run. 

While the MCMC is running, the iteration number, likelihood values, and average 
standard deviation of split frequencies (ASDSF) are printed to the screen. The ASDSF should 
be decreasing toward 0, indicating the tree topologies sampled from the two different runsrafe 
getting similar and converging to the same (stationary) distribution. Typically, ASDSF < 0,01 
is a good sign of consistency between runs. 

To summarize the parameters and trees after the MCMC, use 

sump 

sumt 
By default, the first 25% samples are discarded as burnin. This can be adjusted according to 
the likelihood traces from the two runs (e.g. viewed in Tracer). The effective sample size (ESS) 
is an important indicator to judge if sufficient MCMC samples are collected. Ideally, the ESS 
should be larger than 100 for all parameters. We may need to increase the chain length and 
adjust MCMC proposals to improve the estimates. 

The consensus tree including all fossils is highly unresolved due to the uncertainty in the 
placement of the fossils. In order to display the node ages clearly, we can redraw an extant taxa 
tree by pruning the fossils. The output filename is changed to avoid overwriting the existing ones. 
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delete fossils 

sumt output = hym.rf 
Note here this sumt command does not rerun the MCMC but just for displaying the extant 
taxon tree. 


2.6 Node dating 


In node dating, we calibrate the internal nodes of the tree instead of the fossils above. The 
morphological characters of the fossils are not used, thus we remove them. 
delete fossils 
exclude 24 130 168 
calibrate 
Tenthredinidae = offsetgamma(100,150,25) 
@epSirOruApo = truncatednormal (140,175, 25) 
prset nodeagepr = calibrated 
The birth-death prior under diversified sampling (Hôhna et al., 2011) is used for the time 
tree. Compared to the EBD prior used above, there is no fossil sampling parameter in this case. 


The priors for the root age, speciation and extinction rates are not changed. 
prset brlenspr = clock:bixthdeath 


prset samplestrat = divefsity 


prset sampleprob - 0.0001 
prset speciationpr = exp(10) 
prset extinctionpr = beta(1,1) 
prset treeagepr - offsetexp(300, 390) 
The uniform tree prior (clock: uni form) is also applicable. 


It is important to force the calibrated clade to be monophyletic and enable the.constraints. 

The constraint Hymenoptera helps to root the tree properly. 
prset topologypr = constraint (Hymenoptera, Tenthredimidae, 
CepSirOruApo) 

The output filename is changed to avoid overwriting the existing ones. If the node dating 
analysis is continued after the total-evidence dating in the same MrBayes session, the starting 
values also need to be reset. The other settings are kept the same as in the total-evidence dating 
analysis. 


mcmcp filename = hym.nd startp = reset startt = random 


2.7 Non-clock analysis 


Lastly, we run an analysis without the molecular clock assumption and calibrations. The 
branch lengths are measured by genetic distance (expected substitutions per site). This is a 
typical analysis most people do using MrBayes. We do not use fossils, and do not constrain the 
topology so that they are uniformly distributed. 

delete fossils 


prset brlenspr = uncons:gammadir (1,1,1,1) 
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prset topologypr = uniform 

mcmc filename = hym.un 

sump 

sumt 

The prior for branch lengths is gamma-Dirichlet(1, 1, 1, 1) (Rannala et al., 2012; 

Zhang et al., 2012), which assigns gamma(1, 1) for the tree length and uniform Dirichlet for 
the proportion of branch lengths. The compound Dirichlet prior was shown to help avoid 
overestimating the tree length (Zhang et al., 2012) and is now the default prior in MrBayes 
(since 3.2.3). 


3 Results and discussion 


The posterior estimates of the parameters are summarized in separate files, which can 
be opened using a text editor. The information is also printed to the screen. The partition 
rate multipliers are in hym. *.pstat (here * is to match te and nd, the same below). The 
morphologies (mt fk) and 16S (m(2) ) partitions evolve at similar rate. The 1“ and 2" codon 
positions of Efla (m (3) evolve much more slowly than the 3" codon position (m{ 4}). Thus 
it is reasonable to take into account such significant rate variation across partitions. 

The majority-rule consensus/trees aressummarized in hym. *.con.tre, which can 
be visualized by FigTree (http://tree.bio.ed.acäik/software/figtree) or IcyTree (https://icytree. 
org). The node ages are also in hym. * .vstaf, associated with the bipartition IDs in 
hym.*.parts. The root ID is 0 which includes all taxa, while the ID of Hymenoptera is that 
which excludes the outgroup Raphidioptera ( . 444444444). The extant taxa trees from total- 
evidence dating and node dating under diversified sampling and IGR clock model are shown 
in Fig. 4. The topologies are the same in general, except for the Cephus and Sirex clade. The 
age of Hymenoptera (250 Ma) inferred from total-evidence dating is similar to that from node 
dating, but with a relatively narrower HPD interval. 

The total-evidence dating approach models the fossilization and sampling process 
explicitly, and incorporates different sources of information from the fossil record while 
accounting for the uncertainty of fossil placement. In comparison, the node dating approach 
discards the fossil morphologies, and relies on secondary interpretation of the fossil record as 
node calibrations. Total-evidence dating appears more objective and rigorous, and provides an 
ideal platform for exploring and further improving the models used for Bayesian molecular 
clock dating analysis. 

Comparing the non-clock tree (Fig. 5) with the clock trees (Fig. 4), it is obvious that the 
evolutionary rate is not constant over time. The Xyela and Onycholyda branches evolve much 
more slowly than the Raphidioptera and Orussus/Vespidae clade, and there are indeed dramatic 
rate changes between adjacent branches. In this case, the IGR relaxed clock model is more suitable 
than the autocorrelated TK02 model, and it is not reasonable to assume a strict clock model. 
Apart from these perspectives, a more statistically rigorous model selection procedure (e.g., using 


reversible-jump MCMC or marginal likelihood) needs to be carried out in future studies. 
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A Raphidioptera 
Xyela 
Onycholyda 
Runaria 

Diprion 
Phylacteophaga 
Cephus 

Sirex 

Orussus 


Vespidae 


B Raphidioptera 
Xyela 
Onycholyda 
Runaria 
Diprion 
Phylacteophaga 
Cephus 
Sirex 
Orussus 


Vespidae 


400 350 300 250 200 150 100 50 0 Ma 


Fig. 4 Majority-rule consensus trees of extant taxa from total-evidence dating (A) and node dating (B), 
under the diversified sampling and IGR models 
The node heights are in units of million years and the error bars indicate the 95% HPD,intervals 
The numbers at the internal nodes are the posterior probabilities of the corresponding clades 


Raphidioptera 


Xyela 


Onycholyda 


Runaria 


Diprion 


Phylacteophaga 


Cephus 


Sirex 


Orussus 


Vespidae 


0.1 


Fig. 5 Majority-rule consensus tree of extant taxa from a non-clock analysis under the gamma-Dirichlet prior 
The branch lengths are measured by expected substitutions per site 
The numbers at the internal nodes are the posterior probabilities of the corresponding clades 
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In conclusion, this study provides a brief overview and comparison of total-evidence 
dating and node dating analyses, and demonstrates the functionality of MrBayes using a 
truncated dataset of Hymenoptera. For the analyses and discussion using the full data, please 
see Ronquist et al. (2012a) and Zhang et al. (2016). 


Acknowledgements I sincerely thank Johan Nylander for valuable discussions and for 
organizing a workshop of MrBayes using this tutorial in Stockholm, Sweden. I am grateful to 
Xijun Ni and Rachel C. M. Warnock for insightful review comments. Chi Zhang is supported 
by the 100 Young Talents Program of Chinese Academy of Sciences and partly supported by 
the Strategic Priority Research Program of Chinese Academy of Sciences (XDB26000000). 


MrBaysITHEFZIERF 


Ko W” 
CHERAB AEN sae ARO, PRBS bea ESN DES ISE pă Si ak FL 100044) 
(2 rh [REF EG MIE SPEER JE 100044) 


ME: NAT #|HIMrBayesÙtf†2)-ƒPAZf-WWE2LMItHIEJY. BED FP AN 
TEAS RIE NY ASA H EE, JR SR A Ed EA AEAF, JINNI 
BETTIS WIAA ALLS BAS E AE BOATS 2F, YH PITA RT 
Wr, KONTE TER 

zen: DU RAR AHEM, AP REE, BIRMESE KY 4 MrBayes 


References 


Altekar G, Dwarkadas S, Huelsenbeck J P et al., 2004. Parallel Metropolis coupled Markov chain Monte Carlo for Bayesian 
phylogenetic inference. Bioinformatics, 20: 407-415 

Ayres D L, Darling A, Zwickl D J et al., 2012. BEAGLE: an application programming interface and high-performance 
computing library for statistical phylogenetics. Syst Biol, 61: 170-173 

Drummond A J, Ho S Y W, Phillips M J et al., 2006. Relaxed phylogenetics and dating with confidence. PLoS Biol, 4: 
e88 

Gavryushkina A, Welch D, Stadler T et al., 2014. Bayesian inference of sampled ancestor trees for epidemiology and fossil 
calibration. PLoS Comput Biol, 10: e1003919 

Hastings W K, 1970. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57: 97— 
109 

Heath T A, Huelsenbeck J P, Stadler T, 2014. The fossilized birth-death process for coherent calibration of divergence-time 


estimates. Proc Natl Acad Sci USA, 111: e2957—2966 


201906.00026v1 


chinaXiv 


ChinaXiv 3 F RAT! 


12 


Ho S Y W, Phillips M J, 2009. Accounting for calibration uncertainty in phylogenetic estimation of evolutionary divergence 
times. Syst Biol, 58: 367-380 

Hôhna S, Stadler T, Ronquist F et al., 2011. Inferring speciation and extinction rates under different sampling schemes. Mol 
Biol Evol, 28: 2577-2589 

Huelsenbeck J P, Larget B, Swofford D L, 2000. A compound Poisson process for relaxing the molecular clock. Genetics, 
154: 1879-1892 


Lakner C, van der Mark P, Huelsenbeck J P et al., 2008. Efficiency of Markov chain Monte Carlo tree proposals in Bayesian 
phylogenetics. Syst Biol, 57: 86-103 

Lepage T, Bryant D, Philippe H et al., 2007. A general comparison of relaxed molecular clock models. Mol Biol Evol, 24: 
2669-2680 

Lewis P O, 2001. A likelihood approach to estimating phylogeny from discrete morphological character data. Syst Biol, 50: 
913-925 

Liu E, 2008. BEST: Bayesian estimation of species trees under the coalescent model. Bioinformatics, 24: 2542-2543 

Liu E, Yu, Kubatko L et al., 2009. Coalescent methods for estimating phylogenetic trees. Mol Phylogenet Evol, 53: 320— 
328 

Metropolis N; Rosenbluth A W, Rosenbluth M N et al., 1953. Equation of state calculations by fast computing machines. J 
Chem Phys, 24: 1087-4092 

Rannala B, Zhu T, Yang Z, 202 Tail paradox, partial identifiability, and influential priors in Bayesian branch length 
inference. Mol Biol Evol, 29: 325-335 

Ronquist F, Huelsenbeck J P, 2003. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics, 19: 
1572-1574 

Ronguist F, Klopfstein S, Vilhelmsen L et al., 20124. A4otakevidence approach to dating with fossils, applied to the early 
radiation of the Hymenoptera. Syst Biol, 61: 973-999 

Ronguist F, Teslenko M, van der Mark P et al., 2012b. MrBayes'3.2 efficient Bayesian phylogenetic inference and model 
choice across a large model space. Syst Biol, 61: 539-542 

Stadler T, 2010. Sampling-through-time in birth-death trees. J Theoret Biol, 267: 396-404 

Thorne J L, Kishino H, 2002. Divergence time and evolutionary rate estimation with multilocus data. Syst Biol, 51: 689— 
702 

Van Noorden R, Maher B, Nuzzo R, 2014. The top 100 papers. Nature, 514: 550-553 

Xie W, Lewis PO, Fan Y et al., 2011. Improving marginal likelihood estimation for Bayesian phylogenetic modebselection. 
Syst Biol, 60: 150-160 

Yang Z, 1994a. Estimating the pattern of nucleotide substitution. J Mol Evol, 39: 105-111 

Yang Z, 1994b. Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: 
approximate methods. J Mol Evol, 39: 306-314 

Yang Z, Rannala B, 2006. Bayesian estimation of species divergence times under a molecular clock using multiple fossil 
calibrations with soft bounds. Mol Biol Evol, 23: 212-226 

Zhang C, Rannala B, Yang Z, 2012. Robustness of compound Dirichlet priors for Bayesian inference of branch lengths. Syst 
Biol, 61: 779-784 

Zhang C, Stadler T, Klopfstein S et al., 2016. Total-evidence dating under the fossilized birth-death process. Syst Biol, 65: 
228-249 


ChinaXiv 3 (FRAT! 


Zhang - Molecular clock dating using MrBayes 13 


3 
I 

2 

2 


chinaXiv:201906.00026v1 A 


l| 
m 
inaX 

Chin 


4 Ng. KUA _ ? 
x 


ChinaXiv 3 (FRAT! 


Zhang - Molecular clock dating using MrBayes 15 


3 
I 

2 

2 


chinaXiv:201906.00026v1 A 


1] 
m 
inaX 

Chin 


4 Ng. KUA _ ? 
x 


ChinaXiv 3 (FRAT! 


Zhang - Molecular clock dating using MrBayes 17 


3 
I 

2 

2 


chinaXiv:201906.00026v1 A 


1] 
m 
inaX 

Chin 


4 Ng. KUA _ ? 
x 


ChinaXiv 3 (FRAT) 
853 BU] Wr Y HE zj Wj ®% JR pp. 171-182 


2015441 H VERTEBRATA PALASIATICA figs. 1-3 


È 
N 

2 

2 


chinaXiv:201906.00026v1 A 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


ChinaXiv 3 (FRAT!) 
B536 BU Th A ME zj W ER pp. 183-200 


2015441 H VERTEBRATA PALASIATICA figs. 1-6 


3 
N 

2 

2 


chinaXiv:201906.00026v1 A 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


ChinaXiva (FRAT 
534 SB TH $ HE oh U te TR pp. 201-216 


201541 A VERTEBRATA PALASIATICA figs. 1-2 


3 
N 

2 

2 


chinaXiv:201906.00026v1 A 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


ChinaXiv 3 (FRAT) 
52% 520 T ME oh 1# z ik pp. 217-232 


2014:4H VERTEBRATA PALASIATICA figs. 1-6 


3 
N 

2 

2 


chinaXiv:201906.00026v1 A 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


ChinaXiv 3 ERAT 
?#852% #2 Tr HE zj U 3% JR pp. 233-259 


20144F4J] VERTEBRATA PALASIATICA figs. 1-14 


ka 
y 
o 
2 


chinaXiv:201906.00026v1 A 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


ChinaXiv 3 fF RAF!) 
536 WU ŵ A HE zj W jR pp. 1-264 


20154F1H VERTEBRATA PALASIATICA figs. 1- 


È 
N 

2 

2 


chinaXiv:201906.00026v1 A 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


TỊ 
ivà1FER 
inaXivê 

Chin 


ha Wap. 

& u 

> Ka 
> 


